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(57) ABSTRACT 

A system and a method solve the estimation problem of 
finding reflectance R and illumination L. The system and 
method to solve a functional of the unknown illumination L 
such that a minimum of the functional is assumed to yield a 
good estimate of the illumination L. Having a good estimate 
of the illumination L implies a good estimate of the reflec- 
tance R. The functional uses a variational framework to 
express requirements for the optimal solution. The require- 
ments include: 1) that the illumination Lis spatially smooth; 
2) that the reflectance values are in the interval [0,1] — thus, 
when decomposing the image S, the solution should satisfy 
the constraint L>S; 3) that among all possible solutions, the 
estimate of the illumination L should be as close as possible 
to the image S, so that the contrast of the obtained R is 
maximal; and 4) that the reflectance R complies with typical 
natural image behavior (e.g., the reflectance is piece-wise 
smooth). 
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SYSTEM AND METHOD FOR IMAGE 
ENHANCEMENT, DYNAMIC RANGE 
COMPENSATION AND ILLUMINATION 
CORRECTION 

TECHNICAL FIELD 

[0001] The technical field is enhancement of digital 
images. 

BACKGROUND 

[0002] The human visual system is typically able to adapt 
to lighting variations across scenes, perceiving details in 
regions with very different dynamic ranges. Most image 
recording systems, however, fail this dynamic range com- 
pression task. As a result, images produced by these image 
recording systems are often of poor quality, compared to 
images produced by human perception. Another task that is 
often poorly performed by the image recording systems is 
that of color constancy. Humans perceive color in a way that 
is fairly independent of scene illumination, whereas the 
image recording systems are strongly influenced by spectral 
shifts. 

[0003] The above problems can be stated mathematically 
by describing a relationship between an acquired image S, a 
reflectance of objects of the image R, and an illumination L 
in a pixel -wise multiplication, or: 

S-LR. 

[0004] This expression means that at each point in the 
image S, the color value is the multiplication of the reflec- 
tance value by the illumination value. Given an image S, the 
problem to be solved is removal of the effects of illumination 
and recovery of the reflectance image. That is, given S, find 
both R and L. However, calculation of both R and L is 
typically not possible. The solution then, is to generate a 
methodology that can estimate R and L. 

[0005] Retinex theory deals with compensation for illu- 
mination in images. The goal of the Retinex theory is to 
decompose a given image S into the reflectance image R, 
and the illumination image L, such that at each point (x,y) 
in the image domain, S(x,y) equals the product of R(x,y) and 
L(x,y). The benefits of such decomposition include the 
possibility of removing illumination effects of back/front 
lighting, enhancing photographs and other image capture 
methods that include spatially varying iflurnination, such as 
images that contain indoor and outdoor zones, and correct- 
ing colors in images by removing illumination- induced color 
shifts. 

[0006] As noted above, recovering the illumination L from 
a given image S is known to be a mathematically ill -posed 
problem, and known algorithms vary in the manner and 
effectiveness of overcoming this limitation. The Retinex 
approach provides the framework for one such method. The 
Retinex methodology was motivated by Edward Land's 
research of the human visual system, which is described in 
R. H. Land, "The Retinex Theory of Color Vision/'Sri. 
Amer., Vol. 237, pp. 108-128 (1977). 

[0007] The first Retinex algorithms were of the random 
walk type. Subsequent Retinex algorithms used homomor- 
phic filtering. Yet another group of Retinex algorithms are 
based on solving a Poisson equation. 



[0008] FIG. 1 is a block diagram of an algorithm 10 that 
is representative of the general class of prior- art Retinex 
algorithms. In FIG. 1, an input image S 11 is applied to a 
logarithmetic module 12 to produce a logarithm etic version 
13 of the input image S and its two components, illumination 
L and reflectance R. That is, s=log S, l=log L, and r=log R, 
and thereby, s=l+r. Using the logarithmetic version 13 of the 
input images, an estimator module 14 computes an estimate 
17 of the illumination, designated in FIG. 1 as 1*. The 
estimate 17 (1*) is then combined with the logarithmetic 
version 13 (s) at summer 18 to produce an estimate 19 of the 
reflectance (designated r*). Finally, the estimate 19 (r*) is 
converted from a logarithm to a number (antilogarithm) 
corresponding to the logarithm at expander 20, to produce a 
real number value as an estimate 21 of the reflectance 
(designated as R*). Prior art Retinex algorithms usually 
employ the same process as shown in FIG. 1. 

[0009] The Retinex-based algorithms take several differ- 
ent forms. One such form is the random walk algorithm, 
which is a discrete time random process in which the "next 
pixel position" is chosen randomly from neighbors of the 
current pixel position. Random walk type Retinex algo- 
rithms are variants of the following basic formulation: A 
large number of walkers are initiated at random locations of 
the logarithmetic version 13 (s), adopting a gray-value of 
their initial position. An accumulator image A that has the 
same size as s is initialized to zero. As a walker walks 
around, the walker updates A by adding the value of the 
walker to each position (x,y) that the walker visits. The 
illumination image is obtained by normalizing the accumu- 
lator image A, i.e., the value at each position of the accu- 
mulator image A is divided by the number of walkers that 
visited that position. 

[0010] By using many walkers with long paths, one can 
easily verify that the accumulator value at each position 
converges to a Gaussian average of that position's neigh- 
bors. 

[0011] Another type of Retinex algorithm uses homomor- 
phic filtering, where a low-pass filter is used to reconstruct 
1 from s. Homomorphic Retinex algorithms are based on the 
fact that the reflectance image R corresponds to the sharp 
details in the image (i.e., the edges), whereas the illumina- 
tion Lis expected to be spatially smooth. Then, a reasonable 
guess for 1 is l*-LP{s}, where LP is a convolution with a 
wide Gaussian kernel. Thus, the Retinex algorithm using 
homomorphic filtering actually applies the same process as 
the random walk algorithms by a single direct convolution. 

[0012] Since the illumination L is expected to be spatially 
smooth, the derivative of the illumination L should be close 
to zero. However, considering the assumption that the reflec- 
tance R is piece-wise constant, the derivative of the reflec- 
tance R should vanish almost everywhere, with high values 
along edges of an image. Taking the derivative of the sum 
s=l +r and clipping out high derivative peaks, implies that 
the clipped derivative signal corresponds only to the illu- 
mination L. Poisson equation-type Retinex algorithms that 
rely on the Mondrian world model, use the above assump- 
tions on the reflectance R as a piece-wise constant image. 
Applying the Laplacian, and the following clipping opera- 
tion: 
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[0013] where 

t(Aj)=0 

[0014] otherwise, yields the following Poisson equation 

[0015] A solution to the Poisson equation may involve an 
iterative procedure that effectively inverts the Laplacian 
operator. Improvements to the method involve extracting 
discontinuities from the image gradient magnitude instead of 
the Laplacian to provide better boundary conditions. 

[0016] One solution involves use of an iterative algorithm 
having a "reset*' non-linearity that enforces the constraint 
l^s. The algorithm performs the interactive procedure. 



r M+1 =max{_,— — } 



[0017] where D n is a translation operator, shifting the 
image s by the VI th element of a sequence of spirally decaying 
translation vectors. Removing the max operation yields a 
simplified version 



_ r n + p n [r n ] 
^+1 - ^ 



[0018] The above equation is a simple averaging operation 
that smoothes the images. The non-linear (max) operation 
forces the illumination image 1 to satisfy the constraint l*^s. 

[0019] Despite use of these algorithms, current image 
recoding and image enhancement systems and methods 
cannot produce images that are of sufficient quality to be 
comparable to images as perceived by the human vision 
system. Thus, an improved method and system are required 
to better remove illumination effects and to recover the 
reflectance image. 

SUMMARY 

[0020] A system and a method solve the estimation prob- 
lem of finding reflectance R and illumination L. The system 
and method use a functional of the unknown illumination L 
such that a minimum of the functional is assumed to yield a 
good estimate of the illumination L. Having a good estimate 
of the illumination L implies a good estimate of the reflec- 
tance R. 

[0021] The functional uses a variational framework to 
express requirements for the optimal solution. The require- 
ments include: 1) that the illumination L is spatially smooth; 
2) that the reflectance values are in the interval [0,1] — thus, 
when decomposing the image S, the solution should satisfy 
the constraint L>S; 3) that among all possible solutions, the 
estimate of the illumination L should be as close as possible 
to the image S, so that the contrast of the obtained R is 
maximal; and 4) that the reflectance R complies with typical 
natural image behavior (e.g., the reflectance is piece-wise 
smooth). 



[0022] The four requirements lead to a well-behaved opti- 
mization problem having quadratic programming, which is 
convex, with a single solution. The method and the system 
use a numerical algorithm for the solution, where the algo- 
rithm is based on a multi-resolution decomposition of the 
image S. 

[0023] The method and the system are adaptable to dif- 
ferent image situations. In a first situation, given a color 
image in the RGB color-space, applying the algorithm on the 
three color channels separately produces color correction. 
An alternative embodiment involves converting the input 
image S into a chromatic- luminance color space. The algo- 
rithm is then applied to the luminance alone, resulting in an 
image that preserves the input colors, and improves the local 
contrast 

[0024] A second situation accounts for the fact that 
although the illumination L is assumed to be spatially 
smooth on a global level, the assumption is sometimes 
wrong. This situation is shown dramatically in FIG. 9, 
where the illumination L is expected to have a sharp edge on 
the window boundaries. In situations such as that shown by 
the example of FIG. 9, allowing the illumination L to 
include sharp edges may be enabled by forcing piece-wise 
smoothness on the illumination L, rather than assuming 
global smoothness. 

[0025] A third situation does not use the assumption of 
piece-wise smoothness for the reflectance R. Instead, a 
better quality prior image S is used. The better prior S allows 
use of multi-resolution and uniso tropic considerations. 

[0026] A fourth and final situation accounts for the fact 
that some illumination effect is needed in order to give the 
human user a natural feel of the recorded image. Thus, the 
method and the system may include a process that partially 
returns the illumination such that a final output image 
includes some illumination, but with improved contrast, 

[0027] Using the four requirements enumerated above, the 
variational framework and the algorithm provide a novel 
solution to the problem of estimating illumination. The 
variational framework and the algorithm may be applied in 
many image recording situations, including in digital cam- 
eras as an automatic illumination compensation module, 
aimed at improving image reproduction quality, and in 
scanners and printers, as a special effect that improves the 
visual quality of a scanned/printed image. 

DESCRIPTION OF THE DRAWINGS 

[0028] The detailed description will refer to the following 
drawings, in which like numbers refer to like objects, and in 
which: 

[0029] FIG. 1 is a diagram of a prior art algorithm for 
estimating illumination; 

[0030] FIG. 2 is a block diagram of an improved algo- 
rithm for decomposing an image; 

[0031] FIG. 3 is a diagram of a system that uses the 
improved algorithm of FIG. 2; 

[0032] FIG. 4 is a flowchart showing a process according 
to the algorithm of FIG. 2; 

[0033] FIG. 5 is a flowchart of a sub-routine of the process 
of FIG. 4; 
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[0034] FIG. 6 is a flowchart of a further sub-routine of 
FIG. 4; 

[0035] FIG. 7 is a block diagram of an alternate improved 
algorithm; 

[0036] FIG. 8 is a block diagram of yet another alternate 
improved algorithm; and 

[0037] FIG. 9 is an image triplet showing the effects of 
illumination on an image. 

DETAILED DESCRIPTION 

[0038] The human visual system is typically able to adapt 
to lighting variations across scenes (for example, shadows, 
strong illumination source), perceiving details in regions 
with very different dynamic ranges. Most image recording 
systems (e.g., cameras), however, fail this dynamic range 
compression task. As a result, images produced by these 
image recording systems are often of poor quality, compared 
to images produced by human perception. Another task that 
is often poorly performed by the image recording systems is 
that of color constancy. Humans perceive color in a way that 
is fairly independent of scene illumination, whereas the 
images recording systems are strongly influenced by spectral 
shifts. 

[0039] The above problems can be mathematically formu- 
lated by describing a relationship between an acquired 
image S, a reflectance R of objects of the image, and an 
illumination L in a pixel-wise multiplication, or: 

SmL-K. 

[0040] This expression means that at each point in the 
image S, the color value is the multiplication of the reflec- 
tance value by the illumination value. The effect of this 
relationship can be seen by a simple experiment. A person or 
other object in a room is placed near a window (or other light 
source) during daylight such that a strong light source (i.e., 
bright daylight) illuminates part of an image S to be 
recorded. Thus, the illumination is strong through the win- 
dow and weak inside the room. A photograph is then taken, 
and the recorded image S is decomposed into reflectance R 
and illumination L. The person can barely be seen in the 
image S and the illumination L, but can be seen in the 
reflectance R. 

[0041] Given an image S, the problem to be solved is 
removal of the effects of illumination and recovery of the 
reflectance image. That is, given S, find both R and L. 
However, calculation of either R or L is typically not 
possible. The solution then, is to generate a methodology 
that can estimate R and L. 

[0042] FIG. 2 is a block diagram of an improved algo- 
rithm 100 that may be used to estimate the illumination L 
and the reflectance R. In FIG. 2, an image S 101 is input to 
a logarithm module 102 to produce a logarithm s 103 of the 
image S 101 and its two components, illumination L and 
reflectance R. That is, s-log S, 1-log L, and r-log R, and 
thereby, s-l+r. Using the logarithm s 103, an image pro- 
cessing module 104 uses a Projected Normal Steepest 
Descent or similar algorithm, with multi-resolution process- 
ing, to compute an estimate 107 of the illumination, desig- 
nated in FIG. 2 as I*. The estimate 107 (1*) is then combined 
with the output 103 (s) at summer 108 to produce an 
estimate 109 of the reflectance (designated r*). Finally, the 



estimate 109 (r*) is converted from a logarithm to its 
corresponding base number value at exponential module 
110, to produce a number value as an estimate 111 of the 
reflectance (designated as R*). 

[0043] FIG. 3 is a block diagram of a representative 
component (a digital camera) 112 that uses the system and 
the method described herein. The component 112 includes 
an image capture device 114 that captures an image. The 
image capture device may include lenses, a CCD array, and 
other sub-components. An image converter 115 may be an 
optional subcomponent used for pre-processing of the cap- 
tured image. An image processor 116 may include the 
necessary processing software and algorithms to enhance the 
captured image according to the methods and systems 
described herein. The component 112 may include an image 
memory 117 that stores captured images, both before and 
after processing by the image processor 116. Finally, a 
display device 118, which may be a liquid crystal display 
device, for example, provides a visual display for a user of 
the component 112. 

[0044] To arrive at the solution to the problem of estimat- 
ing R and L, a method and a system for applying the method 
begins by formulating the estimation problem of finding R 
and L from S as an optimization problem. That is, a 
functional of one of the unknowns (either L or R) is provided 
such that a minimum of the functional yields the desired 
result. The formulation uses a variational framework to 
express requirements from the optimal solution. The frame- 
work embodies the following requirements: 

[0045] 1. The illumination L is spatially smooth. 

[0046] 2. The reflectance R is restricted to the unit 
interval, which adds the constraint L^S. Since the 
log function is monotone, l^s. 

[0047] 3. Setting l=Const, where Const is any con- 
stant above the maximal value of s, yields a trivial 
solution that satisfies the two previous requirements. 
Thus, a requirement is added such that the illumina- 
tion image 1 is close to the intensity image s, i.e., the 
value of 1 minimizes a penalty term of the form dist 
(l,s), e.g., the norm (l-s) 2 . 

[0048] 4. The reflectance image s-l-r can be assumed 
to have a high prior probability. One of the simplest 
prior functions used for natural images assigns high 
probability to spatially smooth images. 

[0049] 5. The illumination continues smoothly as a 
constant beyond the image boundaries. This is an 
artificial assumption required for boundary condi- 
tions that would have minor effect on the final 
results. 

[0050] Collecting all the above requirements into one 
expression yields the following penalty functional: 

Minimize: /=Il]-J Q C|V/| 2 +a(/-5) 2 +P|V(/-j)| 2 )(iidfy 

Subject to: fes &nd<V^?>-0 on (1) 

[0051] where Q is the support of the image, 8Q is the 

image boundary, and n is the normal to the boundary, a and 
p are free non-negative real parameters. In the functional 
F[l] f the first penalty term (|^1| 2 ) forces spatial smoothness 
on the illumination image. This choice of smoothness pen- 
alty is natural, keeping in mind that minimizing J(|Vl| 2 )dxdy 
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translates into the Euler-Lagrange (EL) equation Al-0. The 
steepest descent solution to the first penalty term is a 
Gaussian smoothing operation with increasing variance of 
the initial condition. As mentioned in the previous section, 
several authors proposed Gaussian smoothing of s for the 
illumination reconstruction. 

[0052] The second penalty term (l-s) 2 forces a proximity 
between 1 and s. The difference between these images is 
exactly r, which means lhat the norm of r should be small 
(i.e., R tends to white). This term is weighted by the free 
parameter a. The main objective of this term is a regular- 
ization of the problem that makes the second penalty term 
better conditioned. Notice that, in addition, the solution 1 is 
forced to be l^s. 

[0053] The third penalty term is motivated by the Baye- 
sian prior probability theory and encourages the reflectance 
image r to be a "visually pleasing image." In practice, this 
prior probability expression forces r to be spatially smooth. 
The third penalty term is weighted by the free parameter (3. 
More complicated Bayesian expressions may be used allow- 
ing sharp edges, textures, 1/fbehavior, etc. As long as this 
expression is purely quadratic, the above minimization 
problem remains fairly simple. 

[0054] The above-defined problem has a Quadratic Pro- 
gramming (QP) form. The necessary and sufficient condi- 
tions for the minimization problem are obtained with the 
Euler-Lagrange equations: 



V (jc, y) e n 



am 
di 



(2) 



[0060] Using integration by parts, J|VG| 2 — /GAG up to 
the boundary conditions. 

[0061] An alternative approach is the Steepest Descent 
(SD) algorithm, where ^n SD is replaced by a constant value 
//g D , such that: 



[0062] where ^^{A} refers to the greatest eigenvalue of 
the linear operator A. This alternative method saves com- 
putations at the expense of a slightly slower convergence. 

[0063] Finally, projecting onto the constraint l^s is done 
by l J -max(lj, s). In practice, G can be calculated by: 

[0064] where 

G A AA/ M , 

[0065] Similarly, /^sd is given by: 



[0066] where 
AaAJ|V| 2 . 

[0067] The Laplacian may be approximated by a linear 
convolution with the kernel k,* P 



[0055] The differential equation does not have to hold 
when l=s. 

[0056] The minimization problem is QP with respect to the 
unknown image 1. A Projected Normalized Steepest Descent 
(PNSD) algorithm, accelerated by a multi-resolution tech- 
nique, is used to solve this minimization problem. 

[0057] The PNSD algorithm requires the application of a 
Normalized Steepest Descent (NSD) iteration that mini- 
mizes the functional F[l], followed by a projection onto the 
constraints. A NSD iteration has the format: 

W-i-Anbd G 

[0058] where ^ and ^ are the illumination images at step 
j and j-1, respectively, G is the gradient of F[l], and // NSD is 
the optimal line-search step size. For Equation (2), the 
gradient of F[l] is given by: 

[0059] and Mnsd & given by: 



/(a|C| 2 +<H-#|Vq 2 ) 



0 1 0 

1 -4 1 
0 1 0 



[0068] and the integrations are approximated by summa- 
tions: 



fi c i 2 *SZ clfl ' m]2 

n m 



[0069] where G[m, n]=G(mAx, nAy). In order to accom- 
modate the boundary conditions, as given in Equation (1), 
the above convolution is applied on an expanded version of 
the image G. The extension is done by replicating the first 
and last columns and rows. After convolution, the additional 
rows and columns are removed. 

[0070] Although simple, the PNSD algorithm usually con- 
verges slowly. Instead of general acceleration schemes, the 
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fact that the unknown image 1 is assumed to be smooth may 
be used. Specifically, the method applies a multi-resolution 
algorithm that starts by estimating a coarse resolution image 
1, expands the resolution image 1 by interpolation, and uses 
the result as an initialization for the next resolution layer. 
This way, only a few iterations at each resolution are 
required for convergence. 

[0071] Summarizing the above, a block diagram of an 
algorithm 120 for the solution of Equation (1) is shown in 
FIG. 4. The algorithm 120 begins with input block 121, 
where an input to the algorithm 120 is an image s of size [N, 
M], and the two parameters a and p. 

[0072] In initialization block 123, a Gaussian pyramid of 
the image s is computed. The thus-constructed pyramid 
contains p resolution layers (from fine (1) to coarse (p)) with 
the current resolution layer (k) set to p. In block 125, T k 
iterations of the PNSD algorithm are applied to the kth 
resolution layer, until all resolution layers are checked, block 
127. In block 129, the next resolution layer is updated. When 
all resolution layers are processed, the result is the final 
output of the algorithm 120. 

[0073] FIG. 5 is a block diagram showing the initializa- 
tion routine 123 in more detail. In block 131, a counter is 
initialized with 1-1. In block 132, the Gaussian pyramid is 
constructed by smoothing the image with a kernel, such as 
the kernel k PYR : 



_L 1 _L 

16 8 16 



J_ 1 J_ 

16 8 16 



[0074] In block 133, the image is decimated by, for 
example, a 2:1 ratio. The process is repeated (block 134) and 
the counter is incremented (block 135) p times where for this 
section p<l g 2 (min (M, N)). This process produces a 
sequence of images {sj^/ conventionally named the 
"Gaussian Pyramid of s. 1 ' The image Sj is the original image 
s, and Sp is an image with the coarsest resolution for the 
Gaussian pyramid. 

[0075] In block 137, the process sets k=p, i.e., starts at the 
coarsest resolution layer k, and sets the initial condition 
lo«max {sp}. 

[0076] FIG. 6 is a block diagram of the PNSD processing 
main routine 125. The routine 125 starts at block 141, where 
for the ^ resolution layer, the routine 125 calculates: 

[0077] Where: A k -^k^^" 2 (k " 1) ; namely a convolution 
of £ with the kernel k 1JliP as specified above, normalized by 
a factor 2< 2 < k - J » 

[0078] Then, for j=l, . . . ; T k : 

[0079] In block 143, the routine 125 calculates the gradi- 
ent: 



[0080] In block 145, the routine 125 calculates // NS d : 

/i B — <G ( A K G>, 

[0081] Where: 



N M 



[0082] In block 147, the routine 125 completes the NSD 
iteration: 

[0083] In block 149, the result is projected onto the 
constraints: 

[0084] The loop processing solves the intermediate prob- 
lem: 

Minimize: Fi £t />J Rk <|V/| 2 +aa-fi £ ) 2 +p|V(/-5 k )|>^y 
Subject to: l£s x and <Vj;"n*>«0 on da 

[0085] Returning to FIG. 4, in block 127, if k>l, the result 
1^ is up scaled (2:1 ratio) by pixel replication into the new 
1q, the initialization for the next resolution k-1 layer. The 
resolution layer is updated k=k-l, and the algorithm pro- 
ceeds by going again to block 125. If k=l, the result 1^, is 
the final output of the algorithm 120. 

[0086] By setting the parameters a=p=0, and removing 
the constraint l>s, the algorithm 120 reduces to the Homo- 
morphic filtering process that was shown to be equivalent to 
the basic formulation of the random walk algorithms. 

[0087] Using the method and the system described above, 
the solution to Equation (1), and the convergence of the 
numerical algorithm 120 can be shown to be unique. That is, 
the variational optimization problem P, given by 

Minimize: 

[0088] 



Minimize: F t t/1 = J (|V/| 2 +(<*(/ - s k )f +■ y3| V(J - s k )\ l )dxdy 
Subject to:/iJ A and (V7,n) = 0 on dCl 



Subject to: /^s k and <Vl, n >=0 on 3Q 

[0089] with a^0 and, p^0, has a unique solution. 

[0090] The algorithm 120 has thus far been described with 
respect to a single channel. In another embodiment, the 
system and the method are applied to color images. When 
processing color images, one option is to deal with each 
color channel separately. Such channel-by-channel process- 
ing may be referred to as 'RGB Retinex* processing. Treat- 
ing the R, G, and B channels separately usually yields a 
color correction effect. For example, applying RGB Retinex 
processing to a reddish image should modify the fllumina- 
tion in such a way that the red hue is removed and the 
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resulting image is brightened and corrected. Therefore, for 
some image types, RGB Retinex processing improves the 
colors. For other image types, such color correction may 
cause color artifacts that exaggerate color shifts, or reduce 
color saturation. 

[0091] In yet another embodiment, colors are mapped into 
a luminance/chrominance color space, such as HSV. Next, 
the Retinex correction is applied, but only to the luminance 
or V layer, and then the colors are mapped back to the RGB 
domain. This method may be referred to as 'HSV Retinex*. 
Color shifts in such cases are less-likely. The advantage of 
HSV Retinex processing is that only a single channel is 
processed. However, using this embodiment colors are not 
corrected with respect to hue (H), and saturation (S). 

[0092] The reflectance image R obtained by the Retinex 
process is sometimes an over-enhanced image. This can be 
explained by the facts that i) the human visual system 
usually prefers illumination in the image, and that ii) 
removal of all the illumination exposes noise that might 
exist in darker regions of the original image. 

[0093] In still another embodiment, a corrected version of 
the reconstructed illumination is added back to the con- 
structed reflectance image. FIG. 7 describes this operation. 
In FIG. 7, an algorithm 200 computes the illumination 
image L=exp(l) from the intensity image S=exp(s), and the 
reflectance image R=S/L. This is shown where an input 
image S 201 is applied to a Retinex module 202 to produce 
an illumination image L 203. The illumination image L 203 
is tuned by a Gamma correction operation 204 with a free 
parameter y, to obtain a new illumination image L'209. A 
reflectance module 206 computes reflectance R 207 from the 
input image S 201 and the illumination image L 203. The 
new illumination image L'209 is multiplied by the reflec- 
tance R 207 at 210 to produce an output image S*211, where 
S'=L'R. The Gamma correction is performed by 




[0094] where W is the white value (equal to 255 in 8-bit 
images). 

[0095] The final result S' is given, therefore, by: 



S' =L' R=jS 



[0096] For y=1, the whole illumination is added back, and 
therefore S'=S. For y =qo » no illumination is returned, and 
S'=R W, which is the same reflectance image, R, as obtained 
by the algorithm 100, stretched to the interval [0, W]. This 
later case can also be considered as multiplying a maximum 
constant illumination W to the reflectance image R. 

[0097] An analog to adding the illumination to the final 
image can be found in the homomorphic filtering approach. 
A linear filter for the illumination calculation in the log 



domain removes high-pass spatial components of s, yet also 
attenuates the tow-pass components by a factor of Yi (where 
i stands for illumination). This is an analog to a gamma 
correction of the illumination with y^ lt since the equation 
for S' can be written in the form: 



w'\w) Rt 



[0098] and therefore: 



s* -w = -(l-w) + r 
7 

= -(low- pass components) + (high- pass components). 



[0099] In an ideal environment, the illumination L is 
assumed to be piece-wise smooth. That is, L is smooth 
everywhere except at strong discontinuities of the input 
image S. This ideal can be approached in a number of ways. 
First, referring again to FIG. (1), both the illumination L and 
the reflectance R are required to be smooth. However, these 
requirements contradict each other because 1 plus r must 
equal s. Thus 1 and r share the discontinuities. The parameter 
p arbitrates this contradiction and makes sure that appropri- 
ate parts of each discontinuity are apportioned to 1 and r. 

[0100] A second approach ensures that the illumination 1 
smoothness requirement be valid everywhere except over 
discontinuities of s, and that the reflectance smoothness be 
valid over discontinuities of s. 

[0101] This second approach requires computing a weight 
function w(Vs) that obtains values close to 1 over most parts 
of the image s, and close to 0 over strong discontinuities. 
Equation (1) then becomes: 

^^H^s)|7/| a «|/-jR + p(i-H<7j))|V(/- ^ 

[0102] The proj ected steepest descent solution to Equation 
(3) composed of a series of gradient steps is: 

= /; + ii[DfMVs)D£(li)) + Z>>(w(v s^l;)) - a(/; - s) + 



[0103] and projections: 

[0104] Where: 

[0105] ft is a constant 

[0106] D/(,) is a forward derivative in the x-direc- 
tion. 

[0107] D x b Q is a backward derivative in the x-direc- 
tion. 

[0108] D y f (.) is a forward derivative in the y-direc- 
tion. 
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[0109] D y b (.) is a backward derivative in the y-direc- 
tion. 

[0110] FIG. 8 illustrates an algorithm 125a, which is an 
alternative to the algorithm 125 in algorithm 120, that may 
be used to solve the functional (Equation (3)). 

[0111] In block 221, the algorithm 125a calculates w(Vs), 
a k and T k : 

[0112] where Th is a constant threshold, 
[0113] and 

[0114] Note that a k changes for each resolution level-a 
higher o^. corresponds to a course resolution level. 

[0115] In block 223, the algorithm 125a calculates: 

G-DJiw&JfD^U) >Z)/(vv(75)Z)/(lD)-a(l .- j)+ 

*)) 

[0116] In block 225, the algorithm 125a completes the SD 
iteration: 

[0117] In block 227, the result is projected into the con- 
straints 

lj-max {l„J k }. 

[0118] In block 229, the algorithm 125a checks if i=T k , 
and if so, processing according to the algorithm 220 ends. 



In the claims: 

1. An image enhancement method, comprising: 

capturing an image; 

constructing a multi-resolution structure comprising one 
or more resolution layers; 

processing each resolution layer using an iterative algo- 
rithm having at least one iteration; 

projecting each processed resolution layer to a subsequent 
resolution layer; 

up -calling each projected resolution layer to the subse- 
quent resolution layer; and 

using the projected resolution layers to estimate an illu- 
mination image. 

2. The method of claim 1, further comprising, for each of 
one or more iterations: 

calculating a gradient of a penalty functional; and 

computing an optimal line-search step size. 

3. The method of claim 2, wherein the penalty functional 
is given by: 



™ = J ( I V/ 1 2 Hail -s)) 2 +p\VV- \ l )dxdy. 



subject to Us and <V1, n >-0 on dQ; wherein Q is a support 

of the image, 3Q is an image boundary, n is a normal to the 
image boundary, and a and (5 are free non-negative real 
numbers 

4. The method of claim 2, wherein the penalty functional 
is given by: 

^t0=i a( w i (Vj)| V/j 2 +a(/-j)+pw 2 (V5) | V/-V* ftdafy 

where i w 1 and w 2 are non-linear functions of the gradient. 

5. The method of claim 1, wherein the iterative algorithm 
is a Projected Normalized Steepest Descent algorithm. 

6. The method of claim 1, wherein the iterative algorithm 
is a Steepest Descent algorithm. 

7. The method of claim 1, wherein a set of constraints 
comprise that the illumination is greater than the image 
intensity, L>S. 

8. The method of claim 1, further comprising applying 
penalty terms, the penalty terms, comprising: 

that the illumination is spatially smooth; 

that the reflectance is maximized; 

that the reflectance is piece-wise smooth. 

9. The method of claim 1, further comprising: 

computing the reflectance image based on the captured 
image and the estimated illumination image; 

computing a gamma correction factor; 

applying the gamma correction factor to the estimated 
illumination image; and 

multiplying the gamma-corrected illumination image and 
the reflectance image, thereby producing a corrected 
image. 

10. A system, embodied in a computer-readable medium, 
for enhancing digital images, comprising: 

a log module that receives an input digital image S and 
computes a logarithm s of the input digital image; 

an illumination estimator module that produces an esti- 
mate 1* of an illumination component L of the input 
digital image S, wherein the estimator module employs 
a construct comprising one or more resolution layers, 
and an iterative algorithm that processes each of the 
one or more resolution layers; and 

a summing node that sums the logarithm s and a negative 
of the estimate 1* to produce an estimate r* of a 
logarithm of a reflectance component R of the input 
digital image S, wherein a processed resolution layer is 
used to up-scale a subsequent resolution layer. 

11. The system of claim 10, wherein the iterative algo- 
rithm, for each of one or more iterations: 

calculates a gradient of a penalty functional; and 

computes an optimal line-search step size. 

12. The method of claim 11, wherein the penalty func- 
tional is given by: 
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/=■[/]= J ( I V i | 2 Hail - s)) 2 + £ | V (/ - 5) \ 2 )dxdy 



subject to l^s and <V1, n >-0 on dQ; wherein Q is a support 

of the image, dQ is an image boundary, n is a normal to the 
image boundary, and a and p are free non-negative real 
numbers. 

13. The system of claim 10, wherein the penalty func- 
tional is given by: 

^KL-lQ^iCV^Iv/paoaf/^+p^cVj) |v/- 

where Wj and W 2 are non-linear functions of the gradient. 

14. The system of claim 10, wherein the iterative algo- 
rithm is a Projected Normalized Steepest Descent algorithm. 

15. The system of claim 10, wherein the iterative algo- 
rithm is a Steepest Descent algorithm. 

16. The system of claim 10, wherein each of the one or 
more resolution layers is projected onto constraints, and 
wherein the constraints comprise that the illumination is 
greater than the image intensity, L>S; 

17. The system of claim 10, further comprising penalty 
terms, the penalty terms comprising: 

that the illumination is spatially smooth 

that the reflectance is maximized; and 

that the reflectance is piece- wise smooth. 

18. The system of claim 10, further comprising: 

a module that computes reflectance and illumination 
images based on the input digital image and the esti- 
mated illumination image; 



a gamma correction module that computes a gamma 
correction factor and applies the gamma correction 
factor to the estimated illumination image; and 

a node that multiples the gamma-corrected illumination 
image and the reflectance image, thereby producing a 
corrected digital image. 

19. A method for enhancing an image S, the image S 
comprising a reflectance R and an illumination L, the 
method comprising: 

constructing a multi-resolution image structure having 
one or more resolution layers; 

processing the resolution layers using an iterative algo- 
rithm; 

projecting the processed resolution layers onto a set of 
constraints, the set of constraints comprising boundary 
conditions and that L>S; and 

using the projected resolution layers to estimate an illu- 
mination image. 

20. The method of claim 19, wherein the image S is a 
RGB domain color image, the method further comprising: 

mapping colors R, G, B of the image S into a luminance/ 
chrominance color space; 

applying a correction factor to a luminance layer; and 

mapping the luminance/chrominance colors back to the 
RGB domain. 

***** 
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